{
 "cells": [
  {
   "cell_type": "code",
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2025-08-09T11:09:43.495448Z",
     "start_time": "2025-08-09T11:09:43.492692Z"
    }
   },
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import arange\n",
    "\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:09:43.516495Z",
     "start_time": "2025-08-09T11:09:43.513136Z"
    }
   },
   "cell_type": "code",
   "source": [
    "param = np.array([\n",
    "    [30,50,6,18.193],\n",
    "    [30,70,8,9.12],\n",
    "    [30,90,10,9.276],\n",
    "    [40,50,6,30.243],\n",
    "    [40,70,8,23.856],\n",
    "    [40,90,10,35.673],\n",
    "    [50,50,6,6.586],\n",
    "    [50,70,8,6.586],\n",
    "    [50,90,10,8.323],\n",
    "])\n"
   ],
   "id": "df285f9baf5dba76",
   "execution_count": 2,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:09:44.555028Z",
     "start_time": "2025-08-09T11:09:44.552115Z"
    }
   },
   "cell_type": "code",
   "source": "",
   "id": "2b850109c40f3fb3",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T12:42:14.187644Z",
     "start_time": "2025-08-09T12:42:14.170558Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def f(T, H, SC, k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0):\n",
    "\n",
    "    # 变量\n",
    "    solid_content = SC / 100  # 固含量\n",
    "    T = T + 273.15  # 温度，单位为开尔文\n",
    "\n",
    "    # 常量\n",
    "    # 常量\n",
    "    C_0 = 24e-3  # 任意的比例系数，用于确定DMF质量与时间的关系\n",
    "    m_D_0 = 24e-3  # 初始时刻DMF的质量，单位为kg\n",
    "    m_S = 6e-3  # 环丁砜总的质量，单位为g。无论是否溶解，环丁砜不会从体系中消失，总的质量是不变的\n",
    "    m_C = (m_D_0 + m_S) * solid_content  # 初始时刻醋酸纤维素的质量，单位为kg，与初始固含量有关\n",
    "    ro_D = 0.948e3  # DMF的密度，单位为kg/m^3\n",
    "    ro_S = 1.261e3  # 环丁砜的密度，单位为kg/m^3\n",
    "    ro_C = 1.3e3  # 醋酸纤维素的密度，单位为kg/m^3\n",
    "    V = m_D_0 / ro_D + m_S / ro_S + m_C / ro_C  # 总体积，单位为m^3，认为不变\n",
    "    A = 6.09451\n",
    "    B = 2725.96\n",
    "    C = 28.209 # 由温度计算饱和汽压时的三个常数\n",
    "    A_0 = -k_1 * (10 ** (A - B / (T + C)))*133.322 * (1 - k_H * H / 100) / ro_D / V\n",
    "    print(k_1)\n",
    "    print(A)\n",
    "    print(B)\n",
    "    print(C)\n",
    "    print(k_H)\n",
    "    print(ro_D)\n",
    "    print(\"A0\",A_0)\n",
    "    n_correction_factor = 6  # 斯托克斯 - 爱因斯坦方程中的修正系数，对于大分子乃至宏观颗粒是6，对于小分子则体积越小，该值也会越小\n",
    "    viscosity_rate = 0.92e-3  # 在20℃下DMF粘度与水粘度的比值\n",
    "    k = 1.380649e-23  # 玻尔兹曼常数，单位为J/K\n",
    "    r = 0.3413e-9  # 环丁砜的分子半径\n",
    "    N_A = 6.02214076e23  # 阿伏伽德罗常量，单位为mol^-1\n",
    "\n",
    "    # 函数\n",
    "    # DMF在蒸发的过程中，其质量随时间变化的函数\n",
    "    def m_D(t):\n",
    "        return C_0 * np.exp(A_0 * t)\n",
    "\n",
    "\n",
    "    # 未能溶解而析出成液滴的环丁砜的质量随时间变化的函数\n",
    "    def m_S_out(t):\n",
    "        return max(0, m_S - m_D(t) * S_S)\n",
    "\n",
    "\n",
    "    # 未能溶解而析出成固体的醋酸纤维素的质量随时间变化的函数\n",
    "    def m_C_out(t):\n",
    "        return max(0, m_C - m_D(t) * S_C)\n",
    "\n",
    "\n",
    "    # 析出的醋酸纤维素的体积占比随时间变化的函数\n",
    "    def phi_C_out(t):\n",
    "        return m_C_out(t) / ro_C / V\n",
    "\n",
    "\n",
    "    # # 由水在30、40、50（℃）下的粘度和DMF与水粘度的关系计算不同温度下DMF的粘度，单位为mPa*s\n",
    "    # def eta_0(temp):\n",
    "    #     if temp == 30+273.15:\n",
    "    #         return 0.8007 * viscosity_rate\n",
    "    #     elif temp == 40+273.15:\n",
    "    #         return 0.6560 * viscosity_rate\n",
    "    #     elif temp == 50+273.15:\n",
    "    #         return 0.5494 * viscosity_rate\n",
    "\n",
    "\n",
    "    # 不同温度下混合溶液总的粘度随时间变化的函数\n",
    "    def eta(temp, t):\n",
    "        return eta0 * (1 + 2.5 * phi_C_out(t))\n",
    "\n",
    "\n",
    "    # 扩散系数。假定一个小液滴内仅有一个分子，且分子为球形\n",
    "    def D(temp, t):\n",
    "        return k * temp / (n_correction_factor * np.pi * r * eta(temp, t))\n",
    "\n",
    "\n",
    "    # 液滴运动的平均速度\n",
    "    def v(temp, t):\n",
    "        return k_3 * (D(temp, t) ** 0.5)\n",
    "\n",
    "\n",
    "    # 析出的环丁砜的分子数密度随时间变化的函数\n",
    "    def n_molecular_number_density(t):\n",
    "        return m_S_out(t) / (ro_S*4/3*np.pi*r**3) / V\n",
    "\n",
    "\n",
    "    # 小液滴间的碰撞频率\n",
    "    def Z(temp, t):\n",
    "        return (2 ** 0.5) * n_molecular_number_density(t) * np.pi * (r ** 2) * v(temp, t)\n",
    "\n",
    "    tf = 0\n",
    "    for t in np.arange(0, 500, 0.01):\n",
    "\n",
    "        if m_D(t) < 0.024*0.1:\n",
    "            tf = t\n",
    "            break\n",
    "    # print(\"tf：\",tf)\n",
    "\n",
    "    dy = lambda m_s, t: k_4 * Z(T, t) + k_5 * v(T, t) * (m_S_out(t) - m_s) / V * m_s\n",
    "    t = arange(1, tf, 0.01)\n",
    "    sol = odeint(dy, 0, t)\n",
    "\n",
    "    if len(sol.T[0]) == 0:\n",
    "        re = 10000\n",
    "    else:\n",
    "        re = k_6*sol.T[0][-1]\n",
    "    # print(sol.T[0][-1])\n",
    "    return re\n",
    "\n",
    "# print(f(30,50,6,9.3e-9,1e10,8.9e-6,6.15e-6,2.94,3.9e-1,2.31,3.768e-1,99.6))\n",
    "print(f(30.03009834, 75.42591033 , 6.25121123,0.0022, 0.1369, 0.1754, 0.1616, 0.3022, 0.2989, 0.6281, 0.8579, 0.7592))\n",
    "\n"
   ],
   "id": "4598c346b7688cc5",
   "execution_count": 87,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:14:14.728913Z",
     "start_time": "2025-08-09T11:14:14.680578Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def pred(k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0):\n",
    "    predictions = []\n",
    "    for i in range(param.shape[0]):\n",
    "        T = param[i][0]\n",
    "        H = param[i][1]\n",
    "        SC = param[i][2]\n",
    "        re = f(T,H,SC,k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0)\n",
    "        predictions.append(re)\n",
    "    return predictions\n",
    "print(pred(0.01,1,1,1,1,0.9,0.9,0.2,0.9))"
   ],
   "id": "1cf1e6dd5880b7b5",
   "execution_count": 21,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:14:19.338880Z",
     "start_time": "2025-08-09T11:14:19.334602Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# def target(T,H,SC):\n",
    "#     tar = 0\n",
    "#     for i in range(9):\n",
    "#         if param[i,0]==T and param[i,1]==H and param[i,2]==SC:\n",
    "#             tar = param[i,3]\n",
    "#     return tar\n",
    "# print(target(30,50,8))\n",
    "def target():\n",
    "    targets = []\n",
    "    for i in range(param.shape[0]):\n",
    "        t = param[i][-1]\n",
    "        targets.append(t)\n",
    "    return targets\n",
    "print(target())"
   ],
   "id": "7fec16d32cdfff96",
   "execution_count": 22,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:14:27.792979Z",
     "start_time": "2025-08-09T11:14:27.767933Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# def loss(pred, target):\n",
    "#     return np.mean((pred - target)**2)\n",
    "\n",
    "def loss(pred, target):\n",
    "    loss = 0\n",
    "    for i in range(param.shape[0]):\n",
    "        l = np.linalg.norm(pred[i] - target[i])\n",
    "        loss += np.linalg.norm(pred[i] - target[i])\n",
    "    return loss\n",
    "print(loss(pred(0.01646966, 0.03255624, 0.06004486, 0.07621205, 0.97623144,0.85945502, 0.89168254, 0.39898241,0.9),target()))"
   ],
   "id": "cba65906f0368f5a",
   "execution_count": 24,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:14:38.661062Z",
     "start_time": "2025-08-09T11:14:38.657584Z"
    }
   },
   "cell_type": "code",
   "source": [
    "kh_bound = np.array([\n",
    "    [0.7,1],\n",
    "    [0.7,1],\n",
    "    [0,0.8],\n",
    "    [0,1]\n",
    "])\n",
    "\n",
    "k_bound = np.array([\n",
    "    [0.0005,0.03],\n",
    "    [0.001,1],\n",
    "    [0.001,1],\n",
    "    [0.001,1],\n",
    "    [0.001,1],\n",
    "])"
   ],
   "id": "ec0ff9c50f53d774",
   "execution_count": 25,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:40:18.086909Z",
     "start_time": "2025-08-09T11:40:18.083324Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def printre(kkb):\n",
    "    k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0 = kkb\n",
    "    predictions = []\n",
    "    t = target()\n",
    "    for i in range(param.shape[0]):\n",
    "        T = param[i][0]\n",
    "        H = param[i][1]\n",
    "        SC = param[i][2]\n",
    "        re = f(T,H,SC,k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0)\n",
    "        tt = t[i]\n",
    "        print(f\"预测值{re},真实值{tt}\")"
   ],
   "id": "73dd088928c78e1c",
   "execution_count": 47,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:53:57.595085Z",
     "start_time": "2025-08-09T11:53:57.587551Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plot_SA_history():\n",
    "    plt.figure(figsize=(14, 10))\n",
    "\n",
    "    # 子图1：损失函数变化\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.semilogy(history['best_loss'], 'r-', label='Best Loss')\n",
    "    plt.semilogy(history['current_loss'], 'b--', alpha=0.5, label='Current Loss')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Loss')\n",
    "    plt.title('Loss变化曲线')\n",
    "    plt.legend()\n",
    "    plt.grid(True, which=\"both\", ls=\"--\")\n",
    "\n",
    "    # 子图2：温度衰减曲线\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.plot(history['temperature'], 'g-')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('Temperature')\n",
    "    plt.title('温度下降曲线')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图3：k参数演化\n",
    "    plt.subplot(2, 2, 3)\n",
    "    k_array = np.array(history['k_params'])\n",
    "    for i in range(5):\n",
    "        plt.plot(k_array[:, i], label=f'k_{i+1}')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('k参数值')\n",
    "    plt.title('k参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    # 子图4：kh参数演化\n",
    "    plt.subplot(2, 2, 4)\n",
    "    kh_array = np.array(history['kh_params'])\n",
    "    for i in range(4):\n",
    "        plt.plot(kh_array[:, i], label=f'kh_{i+1}')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.ylabel('物理参数值')\n",
    "    plt.title('物理参数变化曲线')\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    plt.grid(True, ls=\"--\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ],
   "id": "b268117688e30fd4",
   "execution_count": 63,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T11:53:58.332049Z",
     "start_time": "2025-08-09T11:53:58.329155Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# 在SA函数外部初始化记录容器\n",
    "history = {\n",
    "    'temperature': [],\n",
    "    'best_loss': [],\n",
    "    'current_loss': [],\n",
    "    'k_params': [],\n",
    "    'kh_params': []\n",
    "}"
   ],
   "id": "73818db199139138",
   "execution_count": 64,
   "outputs": []
  },
  {
   "metadata": {},
   "cell_type": "code",
   "execution_count": null,
   "source": "",
   "id": "1866dbc659900a26",
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T12:04:31.735361Z",
     "start_time": "2025-08-09T12:04:27.434928Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def SA(t0,tf,alpha,iter):\n",
    "    global history\n",
    "    flag = 0\n",
    "    t = t0\n",
    "    k_10,k_30,k_40,k_50,k_60,k_H0,S_C0,S_S0,eta00 = np.array([0.01,0.1,0.1,0.1,0.3,0.3,0.8,0.8,0.8])\n",
    "    k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c = k_10,k_30,k_40,k_50,k_60,k_H0,S_C0,S_S0,eta00\n",
    "\n",
    "    lc = loss(pred(k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c),target())\n",
    "\n",
    "    k_1b,k_3b,k_4b,k_5b,k_6b,k_Hb,S_Cb,S_Sb,eta0b = k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c\n",
    "    lb = 1000\n",
    "    for i in range(iter):\n",
    "        k_1n = k_1c + np.random.normal(0,0.001)\n",
    "        k_1n = np.clip(k_1n,1e-4,0.02)\n",
    "        k_3n = k_3c + np.random.normal(0,0.01)\n",
    "        k_3n = np.clip(k_3n,1e-4,1)\n",
    "        k_4n = k_4c + np.random.normal(0,0.01)\n",
    "        k_4n = np.clip(k_4n,1e-4,1)\n",
    "        k_5n = k_5c + np.random.normal(0,0.01)\n",
    "        k_5n = np.clip(k_5n,1e-4,1)\n",
    "        k_6n = k_6c + np.random.normal(0,0.01)\n",
    "        k_6n = np.clip(k_6n,1e-4,10)\n",
    "        k_Hn = k_Hc + np.random.normal(0,0.01)\n",
    "        k_Hn = np.clip(k_Hn,0.1,0.5)\n",
    "        S_Cn = S_Cc + np.random.normal(0,0.01)\n",
    "        S_Cn = np.clip(S_Cn,0,1)\n",
    "        S_Sn = S_Sc + np.random.normal(0,0.01)\n",
    "        S_Sn = np.clip(S_Sn,0,1)\n",
    "        eta0n = eta0c + np.random.normal(0,0.01)\n",
    "        eta0n = np.clip(eta0n,0,1)\n",
    "\n",
    "        pre = pred(k_1n,k_3n,k_4n,k_5n,k_6n,k_Hn,S_Cn,S_Sn,eta0n)\n",
    "        ln = loss(pre,target())\n",
    "\n",
    "        # for p in pre:\n",
    "        #     if p < 0.001:\n",
    "        #         flag = 1\n",
    "        # if flag==1:\n",
    "        #     continue\n",
    "\n",
    "        if ln < lc or np.random.rand() < np.exp(-(ln-lc)/t):\n",
    "            k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c = k_1n,k_3n,k_4n,k_5n,k_6n,k_Hn,S_Cn,S_Sn,eta0n\n",
    "            lc = ln\n",
    "            if lc < lb:\n",
    "                k_1b,k_3b,k_4b,k_5b,k_6b,k_Hb,S_Cb,S_Sb,eta0b = k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c\n",
    "                lb = lc\n",
    "        flag = 0\n",
    "        kkb = [k_1b,k_3b,k_4b,k_5b,k_6b,k_Hb,S_Cb,S_Sb,eta0b]\n",
    "        t *= alpha\n",
    "        if t < tf:\n",
    "            break\n",
    "        kc = np.array([k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c])\n",
    "        print(f\"iter{i}, lb{lb},kc{kc},ln{ln},\")\n",
    "        k_params = np.array([k_1c,k_3c,k_4c,k_5c,k_6c])\n",
    "        kh_params = np.array([k_Hc,S_Cc,S_Sc,eta0c])\n",
    "        history['temperature'].append(t)\n",
    "        history['best_loss'].append(lb)\n",
    "        history['current_loss'].append(lc)\n",
    "        history['k_params'].append(k_params.copy())\n",
    "        history['kh_params'].append(kh_params.copy())\n",
    "    return kkb\n",
    "re = SA(300,0.001,0.9,1000)\n",
    "print(re)\n",
    "printre(re)\n",
    "plot_SA_history()\n",
    "for key in history.keys():\n",
    "    history[key].clear()"
   ],
   "id": "36d7d95fa58cf227",
   "execution_count": 78,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T09:03:38.798795Z",
     "start_time": "2025-08-09T09:03:38.796237Z"
    }
   },
   "cell_type": "code",
   "source": "",
   "id": "1087da4a08199b1b",
   "execution_count": null,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T10:30:15.537169Z",
     "start_time": "2025-08-09T10:30:15.529062Z"
    }
   },
   "cell_type": "code",
   "source": [
    "# printre([0.01646966, 0.03255624, 0.06004486, 0.07621205, 0.97623144],[0.85945502, 0.89168254, 0.39898241])\n",
    "# 预测值22.987726671540557,真实值18.193\n",
    "# 预测值25.30597950819623,真实值9.12\n",
    "# 预测值28.03341257118123,真实值9.276\n",
    "# 预测值15.312223386779618,真实值30.243\n",
    "# 预测值16.621474421094778,真实值23.856\n",
    "# 预测值18.469303184319546,真实值35.673\n",
    "# 预测值9.989518771611815,真实值6.586\n",
    "# 预测值11.097051579605028,真实值6.586\n",
    "# 预测值12.257962774797742,真实值8.323\n",
    "printre([0.02059353, 0.13302524, 0.16494316, 0.16195733, 0.12423658],[0.8297394 , 0.89689078, 0.28342179])"
   ],
   "id": "13f305ba23439318",
   "execution_count": 281,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T12:04:25.044775Z",
     "start_time": "2025-08-09T12:04:25.041523Z"
    }
   },
   "cell_type": "code",
   "source": [
    "for key in history.keys():\n",
    "    history[key].clear()\n",
    "print(history['temperature'])"
   ],
   "id": "f275ee343e24d313",
   "execution_count": 77,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T12:18:53.400875Z",
     "start_time": "2025-08-09T12:18:53.393639Z"
    }
   },
   "cell_type": "code",
   "source": [
    "list = [0.002177374009341443, 0.13690750443351957, 0.17541473777643418, 0.16163614428118556, 0.30223478158211425, 0.2988828999134941, 0.628070659092061, 0.8578713703472124, 0.759227801695377]\n",
    "list = [round(i,4) for i in list]\n",
    "print(list)"
   ],
   "id": "3d02a510f05f07af",
   "execution_count": 79,
   "outputs": []
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-08-09T10:35:09.735317Z",
     "start_time": "2025-08-09T10:35:09.728658Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "\n",
    "# 真实值和预测值\n",
    "y_true = np.array([18.193, 9.12, 9.276, 30.243, 23.856, 35.673, 6.586, 6.586, 8.323])\n",
    "y_pred = np.array([18.815225, 18.966924, 19.139403, 17.132791, 17.905546, 18.326706, 3.120923, 10.086580, 15.228171])\n",
    "\n",
    "# 计算 MAPE\n",
    "mape = np.mean(np.abs((y_true - y_pred) / y_true)) * 100\n",
    "print(f\"MAPE: {mape:.2f}%\")"
   ],
   "id": "c819be8bea6a19e0",
   "execution_count": 287,
   "outputs": []
  },
  {
   "metadata": {},
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# from scipy.integrate import odeint\n",
    "# import numpy as np\n",
    "# import matplotlib.pyplot as plt\n",
    "# from numpy import arange\n",
    "#\n",
    "# plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体\n",
    "# plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题\n",
    "#\n",
    "# param = np.array([\n",
    "#     [30,50,6,18.193],\n",
    "#     [30,70,8,9.12],\n",
    "#     [30,90,10,9.276],\n",
    "#     [40,50,6,30.243],\n",
    "#     [40,70,8,23.856],\n",
    "#     [40,90,10,35.673],\n",
    "#     [50,50,6,6.586],\n",
    "#     [50,70,8,6.586],\n",
    "#     [50,90,10,8.323],\n",
    "# ])\n",
    "#\n",
    "# def f(T, H, SC, k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0):\n",
    "#\n",
    "#     # 变量\n",
    "#     solid_content = SC / 100  # 固含量\n",
    "#     T = T + 273.15  # 温度，单位为开尔文\n",
    "#\n",
    "#     # 常量\n",
    "#     # 常量\n",
    "#     C_0 = 24e-3  # 任意的比例系数，用于确定DMF质量与时间的关系\n",
    "#     m_D_0 = 24e-3  # 初始时刻DMF的质量，单位为kg\n",
    "#     m_S = 6e-3  # 环丁砜总的质量，单位为g。无论是否溶解，环丁砜不会从体系中消失，总的质量是不变的\n",
    "#     m_C = (m_D_0 + m_S) * solid_content  # 初始时刻醋酸纤维素的质量，单位为kg，与初始固含量有关\n",
    "#     ro_D = 0.948e3  # DMF的密度，单位为kg/m^3\n",
    "#     ro_S = 1.261e3  # 环丁砜的密度，单位为kg/m^3\n",
    "#     ro_C = 1.3e3  # 醋酸纤维素的密度，单位为kg/m^3\n",
    "#     V = m_D_0 / ro_D + m_S / ro_S + m_C / ro_C  # 总体积，单位为m^3，认为不变\n",
    "#     A = 6.09451\n",
    "#     B = 2725.96\n",
    "#     C = 28.209 # 由温度计算饱和汽压时的三个常数\n",
    "#     A_0 = -k_1 * (10 ** (A - B / (T + C)))*133.322 * (1 - k_H * H / 100) / ro_D / V\n",
    "#     n_correction_factor = 6  # 斯托克斯 - 爱因斯坦方程中的修正系数，对于大分子乃至宏观颗粒是6，对于小分子则体积越小，该值也会越小\n",
    "#     viscosity_rate = 0.92e-3  # 在20℃下DMF粘度与水粘度的比值\n",
    "#     k = 1.380649e-23  # 玻尔兹曼常数，单位为J/K\n",
    "#     r = 0.3413e-9  # 环丁砜的分子半径\n",
    "#     N_A = 6.02214076e23  # 阿伏伽德罗常量，单位为mol^-1\n",
    "#\n",
    "#     # 函数\n",
    "#     # DMF在蒸发的过程中，其质量随时间变化的函数\n",
    "#     def m_D(t):\n",
    "#         return C_0 * np.exp(A_0 * t)\n",
    "#\n",
    "#\n",
    "#     # 未能溶解而析出成液滴的环丁砜的质量随时间变化的函数\n",
    "#     def m_S_out(t):\n",
    "#         return max(0, m_S - m_D(t) * S_S)\n",
    "#\n",
    "#\n",
    "#     # 未能溶解而析出成固体的醋酸纤维素的质量随时间变化的函数\n",
    "#     def m_C_out(t):\n",
    "#         return max(0, m_C - m_D(t) * S_C)\n",
    "#\n",
    "#\n",
    "#     # 析出的醋酸纤维素的体积占比随时间变化的函数\n",
    "#     def phi_C_out(t):\n",
    "#         return m_C_out(t) / ro_C / V\n",
    "#\n",
    "#\n",
    "#     # # 由水在30、40、50（℃）下的粘度和DMF与水粘度的关系计算不同温度下DMF的粘度，单位为mPa*s\n",
    "#     # def eta_0(temp):\n",
    "#     #     if temp == 30+273.15:\n",
    "#     #         return 0.8007 * viscosity_rate\n",
    "#     #     elif temp == 40+273.15:\n",
    "#     #         return 0.6560 * viscosity_rate\n",
    "#     #     elif temp == 50+273.15:\n",
    "#     #         return 0.5494 * viscosity_rate\n",
    "#\n",
    "#\n",
    "#     # 不同温度下混合溶液总的粘度随时间变化的函数\n",
    "#     def eta(temp, t):\n",
    "#         return eta0 * (1 + 2.5 * phi_C_out(t))\n",
    "#\n",
    "#\n",
    "#     # 扩散系数。假定一个小液滴内仅有一个分子，且分子为球形\n",
    "#     def D(temp, t):\n",
    "#         return k * temp / (n_correction_factor * np.pi * r * eta(temp, t))\n",
    "#\n",
    "#\n",
    "#     # 液滴运动的平均速度\n",
    "#     def v(temp, t):\n",
    "#         return k_3 * (D(temp, t) ** 0.5)\n",
    "#\n",
    "#\n",
    "#     # 析出的环丁砜的分子数密度随时间变化的函数\n",
    "#     def n_molecular_number_density(t):\n",
    "#         return m_S_out(t) / (ro_S*4/3*np.pi*r**3) / V\n",
    "#\n",
    "#\n",
    "#     # 小液滴间的碰撞频率\n",
    "#     def Z(temp, t):\n",
    "#         return (2 ** 0.5) * n_molecular_number_density(t) * np.pi * (r ** 2) * v(temp, t)\n",
    "#\n",
    "#     tf = 0\n",
    "#     for t in np.arange(0, 500, 0.01):\n",
    "#\n",
    "#         if m_D(t) < 0.024*0.1:\n",
    "#             tf = t\n",
    "#             break\n",
    "#     # print(\"tf：\",tf)\n",
    "#\n",
    "#     dy = lambda m_s, t: k_4 * Z(T, t) + k_5 * v(T, t) * (m_S_out(t) - m_s) / V * m_s\n",
    "#     t = arange(1, tf, 0.01)\n",
    "#     sol = odeint(dy, 0, t)\n",
    "#\n",
    "#     if len(sol.T[0]) == 0:\n",
    "#         re = 10000\n",
    "#     else:\n",
    "#         re = k_6*sol.T[0][-1]\n",
    "#     # print(sol.T[0][-1])\n",
    "#     return re\n",
    "#\n",
    "# # print(f(30,50,6,9.3e-9,1e10,8.9e-6,6.15e-6,2.94,3.9e-1,2.31,3.768e-1,99.6))\n",
    "# print(f(30,50,6,0.01,0.1,0.1,0.1,0.1,0.1,0.9,0.9,0.9))\n",
    "#\n",
    "# def pred(k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0):\n",
    "#     predictions = []\n",
    "#     for i in range(param.shape[0]):\n",
    "#         T = param[i][0]\n",
    "#         H = param[i][1]\n",
    "#         SC = param[i][2]\n",
    "#         re = f(T,H,SC,k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0)\n",
    "#         predictions.append(re)\n",
    "#     return predictions\n",
    "# print(pred(0.01,1,1,1,1,0.9,0.9,0.2,0.9))\n",
    "#\n",
    "# # def target(T,H,SC):\n",
    "# #     tar = 0\n",
    "# #     for i in range(9):\n",
    "# #         if param[i,0]==T and param[i,1]==H and param[i,2]==SC:\n",
    "# #             tar = param[i,3]\n",
    "# #     return tar\n",
    "# # print(target(30,50,8))\n",
    "# def target():\n",
    "#     targets = []\n",
    "#     for i in range(param.shape[0]):\n",
    "#         t = param[i][-1]\n",
    "#         targets.append(t)\n",
    "#     return targets\n",
    "# print(target())\n",
    "#\n",
    "# # def loss(pred, target):\n",
    "# #     return np.mean((pred - target)**2)\n",
    "#\n",
    "# def loss(pred, target):\n",
    "#     loss = 0\n",
    "#     for i in range(param.shape[0]):\n",
    "#         l = np.linalg.norm(pred[i] - target[i])\n",
    "#         loss += np.linalg.norm(pred[i] - target[i])\n",
    "#     return loss\n",
    "# print(loss(pred(0.01646966, 0.03255624, 0.06004486, 0.07621205, 0.97623144,0.85945502, 0.89168254, 0.39898241,0.9),target()))\n",
    "#\n",
    "# def printre(kkb):\n",
    "#     k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0 = kkb\n",
    "#     predictions = []\n",
    "#     t = target()\n",
    "#     for i in range(param.shape[0]):\n",
    "#         T = param[i][0]\n",
    "#         H = param[i][1]\n",
    "#         SC = param[i][2]\n",
    "#         re = f(T,H,SC,k_1,k_3,k_4,k_5,k_6,k_H,S_C,S_S,eta0)\n",
    "#         tt = t[i]\n",
    "#         print(f\"预测值{re},真实值{tt}\")\n",
    "#\n",
    "# def plot_SA_history():\n",
    "#     plt.figure(figsize=(14, 10))\n",
    "#\n",
    "#     # 子图1：损失函数变化\n",
    "#     plt.subplot(2, 2, 1)\n",
    "#     plt.semilogy(history['best_loss'], 'r-', label='Best Loss')\n",
    "#     plt.semilogy(history['current_loss'], 'b--', alpha=0.5, label='Current Loss')\n",
    "#     plt.xlabel('Iteration')\n",
    "#     plt.ylabel('Loss')\n",
    "#     plt.title('Loss变化曲线')\n",
    "#     plt.legend()\n",
    "#     plt.grid(True, which=\"both\", ls=\"--\")\n",
    "#\n",
    "#     # 子图2：温度衰减曲线\n",
    "#     plt.subplot(2, 2, 2)\n",
    "#     plt.plot(history['temperature'], 'g-')\n",
    "#     plt.xlabel('Iteration')\n",
    "#     plt.ylabel('Temperature')\n",
    "#     plt.title('温度下降曲线')\n",
    "#     plt.grid(True, ls=\"--\")\n",
    "#\n",
    "#     # 子图3：k参数演化\n",
    "#     plt.subplot(2, 2, 3)\n",
    "#     k_array = np.array(history['k_params'])\n",
    "#     for i in range(5):\n",
    "#         plt.plot(k_array[:, i], label=f'k_{i+1}')\n",
    "#     plt.xlabel('Iteration')\n",
    "#     plt.ylabel('k参数值')\n",
    "#     plt.title('k参数变化曲线')\n",
    "#     plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "#     plt.grid(True, ls=\"--\")\n",
    "#\n",
    "#     # 子图4：kh参数演化\n",
    "#     plt.subplot(2, 2, 4)\n",
    "#     kh_array = np.array(history['kh_params'])\n",
    "#     for i in range(4):\n",
    "#         plt.plot(kh_array[:, i], label=f'kh_{i+1}')\n",
    "#     plt.xlabel('Iteration')\n",
    "#     plt.ylabel('物理参数值')\n",
    "#     plt.title('物理参数变化曲线')\n",
    "#     plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "#     plt.grid(True, ls=\"--\")\n",
    "#\n",
    "#     plt.tight_layout()\n",
    "#     plt.show()\n",
    "#\n",
    "# # 在SA函数外部初始化记录容器\n",
    "# history = {\n",
    "#     'temperature': [],\n",
    "#     'best_loss': [],\n",
    "#     'current_loss': [],\n",
    "#     'k_params': [],\n",
    "#     'kh_params': []\n",
    "# }\n",
    "#\n",
    "# def SA(t0,tf,alpha,iter):\n",
    "#     global history\n",
    "#     flag = 0\n",
    "#     t = t0\n",
    "#     k_10,k_30,k_40,k_50,k_60,k_H0,S_C0,S_S0,eta00 = np.array([0.01,0.1,0.01,0.1,0.1,0.3,0.9,0.9,0.9])\n",
    "#     k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c = k_10,k_30,k_40,k_50,k_60,k_H0,S_C0,S_S0,eta00\n",
    "#\n",
    "#     lc = loss(pred(k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c),target())\n",
    "#\n",
    "#     k_1b,k_3b,k_4b,k_5b,k_6b,k_Hb,S_Cb,S_Sb,eta0b = k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c\n",
    "#     lb = 1000\n",
    "#     for i in range(iter):\n",
    "#         k_1n = k_1c + np.random.normal(0,0.001)\n",
    "#         k_1n = np.clip(k_1n,1e-4,0.02)\n",
    "#         k_3n = k_3c + np.random.normal(0,0.01)\n",
    "#         k_3n = np.clip(k_3n,1e-4,1)\n",
    "#         k_4n = k_4c + np.random.normal(0,0.01)\n",
    "#         k_4n = np.clip(k_4n,1e-4,1)\n",
    "#         k_5n = k_5c + np.random.normal(0,0.01)\n",
    "#         k_5n = np.clip(k_5n,1e-4,1)\n",
    "#         k_6n = k_6c + np.random.normal(0,0.01)\n",
    "#         k_6n = np.clip(k_6n,1e-4,10)\n",
    "#         k_Hn = k_Hc + np.random.normal(0,0.01)\n",
    "#         k_Hn = np.clip(k_Hn,0.1,0.5)\n",
    "#         S_Cn = S_Cc + np.random.normal(0,0.01)\n",
    "#         S_Cn = np.clip(S_Cn,0,1)\n",
    "#         S_Sn = S_Sc + np.random.normal(0,0.01)\n",
    "#         S_Sn = np.clip(S_Sn,0,1)\n",
    "#         eta0n = eta0c + np.random.normal(0,0.01)\n",
    "#         eta0n = np.clip(eta0n,0,1)\n",
    "#\n",
    "#         pre = pred(k_1n,k_3n,k_4n,k_5n,k_6n,k_Hn,S_Cn,S_Sn,eta0n)\n",
    "#         ln = loss(pre,target())\n",
    "#\n",
    "#         # for p in pre:\n",
    "#         #     if p < 0.001:\n",
    "#         #         flag = 1\n",
    "#         # if flag==1:\n",
    "#         #     continue\n",
    "#\n",
    "#         if ln < lc or np.random.rand() < np.exp(-(ln-lc)/t):\n",
    "#             k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c = k_1n,k_3n,k_4n,k_5n,k_6n,k_Hn,S_Cn,S_Sn,eta0n\n",
    "#             lc = ln\n",
    "#             if lc < lb:\n",
    "#                 k_1b,k_3b,k_4b,k_5b,k_6b,k_Hb,S_Cb,S_Sb,eta0b = k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c\n",
    "#                 lb = lc\n",
    "#         flag = 0\n",
    "#         kkb = [k_1b,k_3b,k_4b,k_5b,k_6b,k_Hb,S_Cb,S_Sb,eta0b]\n",
    "#         t *= alpha\n",
    "#         if t < tf:\n",
    "#             break\n",
    "#         kc = np.array([k_1c,k_3c,k_4c,k_5c,k_6c,k_Hc,S_Cc,S_Sc,eta0c])\n",
    "#         print(f\"iter{i}, lb{lb},kc{kc},ln{ln},\")\n",
    "#         k_params = np.array([k_1c,k_3c,k_4c,k_5c,k_6c])\n",
    "#         kh_params = np.array([k_Hc,S_Cc,S_Sc,eta0c])\n",
    "#         history['temperature'].append(t)\n",
    "#         history['best_loss'].append(lb)\n",
    "#         history['current_loss'].append(lc)\n",
    "#         history['k_params'].append(k_params.copy())\n",
    "#         history['kh_params'].append(kh_params.copy())\n",
    "#     return kkb\n",
    "# re = SA(300,0.001,0.9,1000)\n",
    "# print(re)\n",
    "# printre(re)\n",
    "# plot_SA_history()\n",
    "# for key in history.keys():\n",
    "#     history[key].clear()"
   ],
   "id": "8227df047a1e5ae9",
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
